###########################################################################################
####Replication files for "Experimenting with List Experiments," Public Opinion Quarterly
####Last updated: 3. January 2023
###########################################################################################

################################
################################
################################
####To load all packages necessary for analysis (and more!): 

library(foreign)
library(haven)
library(ggplot2)
library(psych)
library(multilevel)
library(stargazer)
library(xtable)
library(MASS)
library(sandwich)
library(multiwayvcov)
library(prediction)
library(lmtest)
library(tidyverse)
library(data.table)
library(tidyLPA)
library(mclust)
library(plyr)
library(dplyr)
library(knitr)
library(margins)
library(vcd)
library(ggpubr)
library(cowplot)
library(list)
library(emmeans)
library(tidypredict)
library(gmodels)
library(table1)
library(effectsize)
library(MatchIt)
library(RItools)

##Set directory -- set working directory to location of Replication Files folder
setwd("Replication Files")

###To load data:

edweek.df <- read_dta('InPerson_Replication_working.dta')
online.df <- read_dta('Online_Replication_working.dta')

inperson <- subset(edweek.df, white==1) ##This removes respondents who were not white.
inperson <- inperson[!is.na(inperson$list),] ##This removes respondents with missing values on the list experiment.

online <- subset(online.df, white==1) ##This removes respondents who were not white.
online <- online[!is.na(online$list),] ##This removes respondents with missing values on the list experiment.

online_white <- subset(online, int_hisp==0)
online_latinx <- subset(online, int_hisp==1)

inperson_white <- subset(inperson, int_hisp==0)
inperson_latinx <- subset(inperson, int_hisp==1)

inperson <- as.data.frame(inperson)
inperson_white <- as.data.frame(inperson_white)
inperson_latinx <- as.data.frame(inperson_latinx)

online <- as.data.frame(online)
online_white <- as.data.frame(online_white)
online_latinx <- as.data.frame(online_latinx)
online_prime <- subset(online, int_hisp==1 | int_hisp==0)

inperson.latinx <- inperson.white <- inperson
inperson.latinx[, "int_hisp"] <- 1
inperson.white[, "int_hisp"] <- 0

online.latinx <- online.white <- online
online.latinx[, "int_hisp"] <- 1
online.white[, "int_hisp"] <- 0

women_inperson <- inperson[inperson$female == 1,]
men_inperson <- inperson[inperson$female == 0,]

women_online <- online[online$female == 1,]
men_online <- online[online$female == 0,]

women_inperson.latinx <- women_inperson.white <- women_inperson
women_inperson.latinx[, "int_hisp"] <- 1
women_inperson.white[, "int_hisp"] <- 0

men_inperson.latinx <- men_inperson.white <- men_inperson
men_inperson.latinx[, "int_hisp"] <- 1
men_inperson.white[, "int_hisp"] <- 0


################################################################################
##Main Tables and Figures
##Results in Article Text
################################################################################

##Prepare variables
set.seed(1)
describe(inperson$list)
describe(inperson$treat)

inperson$list <- as.integer(inperson$list)
inperson$treat <- as.integer(inperson$treat)


##Table 1. Observed Data from Face-to-Face and Online List Experiments
inperson_table <- table(inperson$list, inperson$treat)
inperson_table
prop.table(inperson_table, 2)

CrossTable(inperson$list, inperson$treat, prop.c=TRUE, prop.r=FALSE, prop.t=FALSE, prop.chisq=FALSE)


online_table <- table(online$list, online$treat)
online_table
prop.table(online_table, 2)

CrossTable(online$list, online$treat, prop.c=TRUE, prop.r=FALSE, prop.t=FALSE, prop.chisq=FALSE)

table(online$int_hisp) 
table(online$int_fem)
table(online$int_control)
table(online$treat)

##Table 2. Difference-in-Differences Results: Face-to-Face Study
reg_result_edweek <- lm(list ~ treat*int_hisp + age + college_grad + female + int_fem, data = inperson)
summary(reg_result_edweek)
stargazer(reg_result_edweek)

emm_edweek = emmeans(reg_result_edweek, specs=revpairwise ~ treat:int_hisp)
emm_edweek$emmeans
emm_edweek$contrasts %>% summary(infer=TRUE) %>% as.data.frame()

##Table 3. Difference-in-Differences by Participant Gender: Face-to-Face Study
reg_result_edweek_women <- lm(list ~ treat*int_hisp + age + college_grad + int_fem, data = women_inperson)
summary(reg_result_edweek_women)

emm_edweek_women = emmeans(reg_result_edweek_women, specs=revpairwise ~ treat:int_hisp)
emm_edweek_women$emmeans
emm_edweek_women$contrasts %>% summary(infer=TRUE) %>% as.data.frame()

reg_result_edweek_men <- lm(list ~ treat*int_hisp + age + college_grad + int_fem, data = men_inperson)
summary(reg_result_edweek_men)

emm_edweek_men = emmeans(reg_result_edweek_men, specs=revpairwise ~ treat:int_hisp)
emm_edweek_men$emmeans
emm_edweek_men$contrasts %>% summary(infer=TRUE) %>% as.data.frame()

##Table 4. Difference-in-Differences Results: Online Study
reg_result_online <- lm(list ~ treat*int_hisp + age + college_grad + female + int_fem, data = online)
summary(reg_result_online)
stargazer(reg_result_online)

emm_online = emmeans(reg_result_online, specs=revpairwise ~ treat:int_hisp)
emm_online$emmeans
emm_online$contrasts %>% summary(infer=TRUE) %>% as.data.frame()

##Table 5. Difference-in-Differences Results by Participant Gender: Online Study
#Online result by gender of respondent
reg_result_online_women <- lm(list ~ treat*int_hisp + age + college_grad + int_fem, data = women_online)
summary(reg_result_online_women)

reg_result_online_women = emmeans(reg_result_online_women, specs=revpairwise ~ treat:int_hisp)
reg_result_online_women$emmeans
reg_result_online_women$contrasts %>% summary(infer=TRUE) %>% as.data.frame()

reg_result_online_men <- lm(list ~ treat*int_hisp + age + college_grad + int_fem, data = men_online)
summary(reg_result_online_men)

reg_result_online_men = emmeans(reg_result_online_men, specs=revpairwise ~ treat:int_hisp)
reg_result_online_men$emmeans
reg_result_online_men$contrasts %>% summary(infer=TRUE) %>% as.data.frame()

##Figure 1. Results Are Robust to Model Type.
diff.in.means.results_edweek <- ictreg(list ~ 1, data = inperson,
                                       treat = "treat", J=4, method = "lm")
summary(diff.in.means.results_edweek)


diff.in.means.results_white <- ictreg(list ~ 1, data = inperson_white,
                                      treat = "treat", J=4, method = "lm")
diff.in.means.results_latinx <- ictreg(list ~ 1, data = inperson_latinx,
                                       treat = "treat", J=4, method = "lm")
summary(diff.in.means.results_white)
summary(diff.in.means.results_latinx)

lm.results <- ictreg(list ~ age + college_grad + female + int_fem + int_hisp, data = inperson,
                     treat = "treat", J=4, method = "lm")
summary(lm.results)

avg.pred.white.lm <- predict(lm.results, newdata = inperson.white, se.fit = TRUE, avg = TRUE) 
avg.pred.latinx.lm <- predict(lm.results, newdata = inperson.latinx, se.fit = TRUE, avg = TRUE)
#avg.pred.diff.lm <- predict(lm.results, newdata = inperson.white, newdata.diff = inperson.latinx, se.fit = TRUE, avg = TRUE, interval="confidence")
avg.pred.diff.lm <- avg.pred.white.lm$fit-avg.pred.latinx.lm$fit

lm.results_white <- ictreg(list ~ age + college_grad + female + int_fem, data = inperson_white,
                           treat = "treat", J=4, method = "lm")
lm.results_latinx <- ictreg(list ~ age + college_grad + female + int_fem, data = inperson_latinx,
                            treat = "treat", J=4, method = "lm")
summary(lm.results_white)
summary(lm.results_latinx)

avg.pred.white.lm <- predict(lm.results, newdata = inperson.white, se.fit = TRUE, avg = TRUE, interval="confidence")
avg.pred.latinx.lm <- predict(lm.results, newdata = inperson.latinx, se.fit = TRUE, avg = TRUE, interval="confidence")
plot(c(avg.pred.white.lm, avg.pred.latinx.lm), labels = c("White Interviewers", "Latinx Interviewers"), main="LM Results")

ml.results <- ictreg(list ~ age + college_grad + female + int_fem + int_hisp, data = inperson, treat = "treat", J=4, method = "ml")
summary(ml.results)

avg.pred.white.ml <- predict(ml.results, newdata = inperson.white, se.fit = TRUE, avg = TRUE)
avg.pred.latinx.ml <- predict(ml.results, newdata = inperson.latinx, se.fit = TRUE, avg = TRUE)
avg.pred.diff.ml <- predict(ml.results, newdata = inperson.white, newdata.diff = inperson.latinx, se.fit = TRUE, avg = TRUE, interval = "confidence")
plot(avg.pred.diff.ml, labels = c("White Interviewers", "Latinx Interviewers", "Difference"), main="ML Results")


ml.ceilingfloor.results <- ictreg(list ~ age + college_grad + female + int_fem + int_hisp, data = inperson, treat = "treat", J=4, method = "ml",
                                  floor=TRUE, ceiling=TRUE, ceiling.fit = "bayesglm",ceiling.formula = ~ age + college_grad + female + int_fem + int_hisp, floor.fit = "bayesglm", floor.formula = ~ age + college_grad + female + int_fem + int_hisp)
summary(ml.ceilingfloor.results)
#avg.pred.white.mle <- predict(ml.constrained.results,
#                              newdata = inperson.white, avg = TRUE, interval = "confidence")
#avg.pred.latinx.mle <- predict(ml.constrained.results,
#                               newdata = inperson.latinx, avg = TRUE, interval = "confidence")
avg.pred.diff.mle <- predict(ml.ceilingfloor.results,
                             newdata = inperson.white, newdata.diff = inperson.latinx,
                             se.fit = TRUE, avg = TRUE, interval="confidence")
plot(avg.pred.diff.mle, labels = c("White Interviewer", "Latinx Interviewer", "Difference"), main="ML Results with Ceiling and Floor Adjustments")




nls.results <- ictreg(list ~ age + college_grad + female + int_fem + int_hisp, data = inperson,
                      treat = "treat", J=4, method = "nls")
summary(nls.results)

avg.pred.white.nls <- predict(nls.results, newdata = inperson.white, se.fit = TRUE, avg = TRUE)
avg.pred.latinx.nls <- predict(nls.results, newdata = inperson.latinx, se.fit = TRUE, avg = TRUE)
avg.pred.diff.nls <- predict(nls.results, newdata = inperson.white, newdata.diff = inperson.latinx, se.fit = TRUE, avg = TRUE, interval = "confidence")
plot(avg.pred.diff.nls, labels = c("White Interviewer", "Latinx Interviewer", "Difference"))


plot(c(avg.pred.diff.ml, avg.pred.diff.nls, avg.pred.diff.mle), labels = c("ML Results", "NLS Results", "Floor/Ceiling"))

fit_ml <- avg.pred.diff.ml$fit
se_ml <- avg.pred.diff.ml$se.fit

interview <- c("White Interviewer", "Latinx Interviewer", "Difference")
model_type <- c("ML", "ML", "ML")

fit_ml <- cbind(fit_ml, se_ml, interview, model_type)
names(fit_ml)[names(fit_ml) == "se_ml"] <- "se"

fit_nls <- avg.pred.diff.nls$fit
se_nls <- avg.pred.diff.nls$se.fit

model_type <- c("NLS", "NLS", "NLS")

fit_nls <- cbind(fit_nls, se_nls, interview, model_type)
names(fit_nls)[names(fit_nls) == "se_nls"] <- "se"

fit_mle <- avg.pred.diff.mle$fit
se_mle <- avg.pred.diff.mle$se.fit

model_type <- c("ML with Floor/Ceiling", "ML with Floor/Ceiling", "ML with Floor/Ceiling")

fit_mle <- cbind(fit_mle, se_mle, interview, model_type)
names(fit_mle)[names(fit_mle) == "se_mle"] <- "se"

predplots_diff.df <- rbind(fit_ml, fit_nls, fit_mle)
predplots_diff.df$model_type <- factor(predplots_diff.df$model_type, levels=c("ML", "NLS", "ML with Floor/Ceiling"))

ggplot(data = predplots_diff.df, aes(x=interview, y=fit, ymin=lwr, ymax=upr, colour=model_type, group=model_type)) + 
  expand_limits(y=c(0,1)) + 
  geom_point(size=8, stat="identity", position=position_dodge(0.4)) +
  geom_errorbar(aes(), width=.1, position=position_dodge(0.4)) +
  scale_colour_manual(values=c("#2c7bb6",  "#5e3c99", "#d7191c"), name="Model Type") + 
  scale_x_discrete(name="\nExperimental Condition",
                   limits=c("White Interviewer", "Latinx Interviewer", "Difference"))+  
  ylab("Estimated Proportion") +
  geom_text(aes(label=format(fit, digits=2)), fontface="bold", colour="white", size=3, position=position_dodge(0.4))+
  geom_hline(yintercept=0, linetype="dashed", color = "red") +
  theme_light(base_size=14) +
  theme(legend.position = c(.9, .75), legend.justification = c(.8, .2))


##Figure 2 is a screen shot from the online study.


##Figure 3. Effects of Latinx Interviewers on Direct Political Attitudes.
##Figures for explicit question
##Face-to-face
##In general, illegal immigrants are more prone to violence
reg_explicit1_inperson <- lm(Q11_1 ~ int_hisp + age + college_grad + female + int_fem, data = inperson)
summary(reg_explicit1_inperson)

##Given the current immigration situation, denying some basic constitutional rights is justified
reg_explicit3_inperson <- lm(Q11_3 ~ int_hisp + age + college_grad + female + int_fem, data = inperson)
summary(reg_explicit3_inperson)

##Online
##In general, illegal immigrants are more prone to violence
reg_explicit1_online <- lm(Q18_1 ~ int_hisp + age + college_grad + female + int_fem, data = online)
summary(reg_explicit1_online)

##Given the current immigration situation, denying some basic constitutional rights is justified
reg_explicit3_online <- lm(Q18_3 ~ int_hisp + age + college_grad + female + int_fem, data = online)
summary(reg_explicit3_online)

modelcoef1 <- tidy(reg_explicit1_inperson) %>% filter(term == "int_hisp") %>% mutate(source = "Face-to-Face", var = "1")
modelcoef2 <- tidy(reg_explicit3_inperson) %>% filter(term == "int_hisp") %>% mutate(source = "Face-to-Face", var = "2")
modelcoef3 <- tidy(reg_explicit1_online) %>% filter(term == "int_hisp") %>% mutate(source = "Online", var = "1")
modelcoef4 <- tidy(reg_explicit3_online) %>% filter(term == "int_hisp") %>% mutate(source = "Online", var = "2")

explicit_fig.df <- rbind(modelcoef1, modelcoef2, modelcoef3, modelcoef4)

labels <- c("1"="Illegal immigrants more \nprone to violence", "2"="Deny illegal immigrants \nbasic rights", "3"="Voted for Trump")

ggplot(data = explicit_fig.df, aes(x=var, y=estimate, ymin=estimate-std.error*qnorm(.975), ymax=estimate+std.error*qnorm(.975), group=source, shape=source)) + 
  geom_hline(yintercept=0, linetype="dashed", colour="red")+
  geom_point(size=4, stat="identity", position=position_dodge(0.2)) +
  geom_errorbar(aes(), width=.1, position=position_dodge(.2)) +
  #geom_text(aes(label=format(coef_result, digits=2)), fontface="bold", colour="white", size=1.8, position=position_dodge(0.2))+
  scale_x_discrete(name="", labels=labels)+
  ylab("Effect of Latinx Interviewer") +
  ggtitle("Coefficient Plot") + 
  coord_flip()+
  theme_bw(base_size = 14) +
  theme(legend.title = element_blank())


################################################################################
##Appendix Results
################################################################################

##Table A1. Face-to-Face Study: Response Rates by Interviewer Condition
##Response rates by interviewer ethnicity
describeBy(edweek.df$agree_response,list(edweek.df$int_hisp, edweek.df$int_fem), mat=TRUE, digits=2) #two grouping variables
describeBy(edweek.df$agree_response, list(edweek.df$int_hisp), mat=TRUE, digits=2)
describeBy(edweek.df$agree_response, list(edweek.df$int_fem), mat=TRUE, digits=2)

##Other calculations
describe(edweek.df$agree_response)

t.test(agree_response ~ int_fem, data=edweek.df)
t.test(agree_response ~ int_hisp, data=edweek.df)
t.test(agree_response ~ int_fem, data=edweek.df, subset=int_hisp==0)
t.test(agree_response ~ int_fem, data=edweek.df, subset=int_hisp==1)
t.test(agree_response ~ int_hisp, data=edweek.df, subset=int_fem==0)

table(inperson$int_hisp) 
table(inperson$int_fem)
table(inperson$treat)

##Descriptive Statistics
##Table A2
stargazer(inperson[c("age", "college_grad", "female", "pid7", "treat", "int_hisp", "int_fem")],
          #type="latex", out="Tables/yougov_summary.tex", float = FALSE, #style="apsr",
          covariate.labels = c("Age", "College Graduate", "Female", "Party Identification", "Immigration Treatment Condition", "Latinx Interviewer", "Female Interviewer") 
)

##Table A3
stargazer(online[c("age", "college_grad", "female", "pid7", "treat", "int_control", "int_hisp", "int_fem")],
          #type="latex", out="Tables/yougov_summary.tex", float = FALSE, #style="apsr",
          covariate.labels = c("Age", "College Graduate", "Female", "Party Identification", "Immigration Treatment Condition", "No Prime Control", "Latinx Prime", "Female Prime") 
)


##Table A4. Omnibus Balance Tests

##Table A4, Row 1
xBalance(int_hisp~ white + pid7 + educ + female + age,
         data=inperson,
         report=c("chisquare.test", "std.diffs"))

##Table A4, Row 2
xBalance(int_fem~ white + pid7 + educ + female + age,
         data=inperson,
         report=c("chisquare.test", "std.diffs"))

##Table A4, Row 3
xBalance(int_hisp~ white + pid7 + educ + female + age,
         data=online_prime,
         report=c("chisquare.test", "std.diffs"))

##Table A4, Row 4
xBalance(int_fem~ white + pid7 + educ + female + age,
         data=online_prime,
         report=c("chisquare.test", "std.diffs"))



##Design Effect Tests
##Table A5. Test for List Experiment Design Effects: Face-to-Face Study
test.value.edweek <- ict.test(inperson$list, inperson$treat, J = 4, gms = TRUE, pi.table = TRUE)
print(test.value.edweek)

test.value.edweek_white <- ict.test(inperson_white$list, inperson_white$treat, J = 4, gms = TRUE, pi.table = TRUE)
print(test.value.edweek_white)

test.value.edweek_latinx <- ict.test(inperson_latinx$list, inperson_latinx$treat, J = 4, gms = TRUE, pi.table = TRUE)
print(test.value.edweek_latinx)

##Table A6. Test ofr List Experiment Design Effects: Online Survey
test.value.online <- ict.test(online$list, online$treat, J = 4, gms = TRUE, pi.table = TRUE)
print(test.value.online)


##Table A7. Face-to-Face Study: Linear Model
reg_result_edweek <- lm(list ~ treat*int_hisp + age + college_grad + female + int_fem, data = inperson)
summary(reg_result_edweek)
stargazer(reg_result_edweek)

##Table A8. Face-to-Face Study: Comparison of Models
lm.results <- ictreg(list ~ age + college_grad + female + int_fem + int_hisp, data = inperson,
                     treat = "treat", J=4, method = "lm")
summary(lm.results)

ml.results <- ictreg(list ~ age + college_grad + female + int_fem + int_hisp, data = inperson, treat = "treat", J=4, method = "ml")
summary(ml.results)

nls.results <- ictreg(list ~ age + college_grad + female + int_fem + int_hisp, data = inperson,
                      treat = "treat", J=4, method = "nls")
summary(nls.results)

ml.ceilingfloor.results <- ictreg(list ~ age + college_grad + female + int_fem + int_hisp, data = inperson, treat = "treat", J=4, method = "ml",
                                  floor=TRUE, ceiling=TRUE, ceiling.fit = "bayesglm",ceiling.formula = ~ age + college_grad + female + int_fem + int_hisp, floor.fit = "bayesglm", floor.formula = ~ age + college_grad + female + int_fem + int_hisp)
summary(ml.ceilingfloor.results)


##Table A9
reg_result_edweek_women <- lm(list ~ treat*int_hisp + age + college_grad + int_fem, data = women_inperson)
summary(reg_result_edweek_women)

emm_edweek_women = emmeans(reg_result_edweek_women, specs=revpairwise ~ treat:int_hisp)
emm_edweek_women$emmeans
emm_edweek_women$contrasts %>% summary(infer=TRUE) %>% as.data.frame()

reg_result_edweek_men <- lm(list ~ treat*int_hisp + age + college_grad + int_fem, data = men_inperson)
summary(reg_result_edweek_men)

emm_edweek_men = emmeans(reg_result_edweek_men, specs=revpairwise ~ treat:int_hisp)
emm_edweek_men$emmeans
emm_edweek_men$contrasts %>% summary(infer=TRUE) %>% as.data.frame()

reg_result_edweek_genderinteract <- lm(list ~ treat*int_hisp*female + age + college_grad + int_fem, data = inperson)
summary(reg_result_edweek_genderinteract)

stargazer(reg_result_edweek_men, reg_result_edweek_women)

vars.order <- c("treat", "int_hisp", "treat:int_hisp", "age", "college_grad", "int_fem", "female", "treat:female", "int_hisp:female", "treat:int_hisp:female")
stargazer(reg_result_edweek_men, reg_result_edweek_women, reg_result_edweek_genderinteract,
          order=paste0("^", vars.order , "$"),
          column.labels = c("\\thead{Male Participants}", "\\thead{Female Participants}", "\\thead{Triple Interaction}"),
          #type="latex", out="Tables/yougov_summary.tex", float = FALSE, #style="apsr",
          covariate.labels = c("Immigration Condition", "Latinx Interviewer", "Immigration Condition*Latinx Interviewer", "Age", "College Graduate", "Female Interviewer", "Female Participant", "Immigration Condition*Female Participant", "Latinx Interviewer*Female Participant", "Immigration Condition*Latinx Interviewr*Female Participant") 
)


##Table A10. Face-to-Face Study: Linear Model with Interviewer Gender Interaction
reg_result_edweek_intfem <- lm(list ~ treat*int_fem + age + college_grad + female + int_hisp, data = inperson)
summary(reg_result_edweek_intfem)
stargazer(reg_result_edweek_intfem)

##Table A11. Difference-in-Differences Results by Interviwer Gender: Face-to-Face Study
emm_edweek_intfem = emmeans(reg_result_edweek_intfem, specs=revpairwise ~ treat:int_fem)
emm_edweek_intfem$emmeans
emm_edweek_intfem$contrasts %>% summary(infer=TRUE) %>% as.data.frame()


##Table A12. Online Study: Linear Model
reg_result_online <- lm(list ~ treat*int_hisp + age + college_grad + female + int_fem, data = online)
summary(reg_result_online)
stargazer(reg_result_online)

emm_online = emmeans(reg_result_online, specs=revpairwise ~ treat:int_hisp)
emm_online$emmeans
emm_online$contrasts %>% summary(infer=TRUE) %>% as.data.frame()

##Table A13. Online Survey: No Prime
reg_result_online_control <- lm(list ~ treat + age + college_grad + female, data = subset(online, int_control==1))
summary(reg_result_online_control)
stargazer(reg_result_online_control)

##Table A14: Online Survey: Estimated Values in No Prime Condition
emm_online_control = emmeans(reg_result_online_control, specs=revpairwise ~ treat, data =subset(online, int_control==1) )
emm_online_control$emmeans
emm_online_control$contrasts %>% summary(infer=TRUE) %>% as.data.frame()


##Table A15. Online Study: Comparison of Models
##Online
lm.results_online <- ictreg(list ~ age + college_grad + female + int_fem + int_hisp, data = online,
                            treat = "treat", J=4, method = "lm")
summary(lm.results_online)

ml.results_online <- ictreg(list ~ age + college_grad + female + int_fem + int_hisp, data = online,
                            treat = "treat", J=4, method = "ml")
summary(ml.results_online)

nls.results_online <- ictreg(list ~ age + college_grad + female + int_fem + int_hisp, data = online,
                             treat = "treat", J=4, method = "nls")
summary(nls.results_online)

ml.ceilingfloor.results_online <- ictreg(list ~ age + college_grad + female + int_fem + int_hisp, data = online, treat = "treat", J=4, method = "ml",
                                         floor=TRUE, ceiling=TRUE, ceiling.fit = "bayesglm",ceiling.formula = ~ age + college_grad + female + int_fem + int_hisp, floor.fit = "bayesglm", floor.formula = ~ age + college_grad + female + int_fem + int_hisp)
summary(ml.ceilingfloor.results_online)


##Figure A1. Online Survey: Results by Model Type
avg.pred.white.ml_online <- predict(ml.results_online, newdata = online.white, se.fit = TRUE, avg = TRUE)
avg.pred.latinx.ml_online <- predict(ml.results_online, newdata = online.latinx, se.fit = TRUE, avg = TRUE)
avg.pred.diff.ml_online <- predict(ml.results_online, newdata = online.white, newdata.diff = online.latinx, se.fit = TRUE, avg = TRUE, interval = "confidence")
plot(avg.pred.diff.ml_online, labels = c("White Interviewer", "Latinx Interviewer", "Difference"), main="Online Study")


#avg.pred.white.mle <- predict(ml.constrained.results,
#                              newdata = inperson.white, avg = TRUE, interval = "confidence")
#avg.pred.latinx.mle <- predict(ml.constrained.results,
#                               newdata = inperson.latinx, avg = TRUE, interval = "confidence")
avg.pred.diff.mle_online <- predict(ml.ceilingfloor.results_online,
                                    newdata = online.white, newdata.diff = online.latinx,
                                    se.fit = TRUE, avg = TRUE, interval="confidence")
plot(avg.pred.diff.mle_online, labels = c("White Interviewer", "Latinx Interviewer", "Difference"), main="ML Results with Ceiling and Floor Adjustments")



avg.pred.white.nls_online <- predict(nls.results_online, newdata = online.white, se.fit = TRUE, avg = TRUE)
avg.pred.latinx.nls_online <- predict(nls.results_online, newdata = online.latinx, se.fit = TRUE, avg = TRUE)
avg.pred.diff.nls_online <- predict(nls.results_online, newdata = online.white, newdata.diff = online.latinx, se.fit = TRUE, avg = TRUE, interval = "confidence")
plot(avg.pred.diff.nls_online, labels = c("White Interviewer", "Latinx Interviewer", "Difference"))


fit_ml_online <- avg.pred.diff.ml_online$fit
se_ml_online <- avg.pred.diff.ml_online$se.fit

interview <- c("White Interviewer", "Latinx Interviewer", "Difference")
model_type <- c("ML", "ML", "ML")

fit_ml_online <- cbind(fit_ml_online, se_ml_online, interview, model_type)
names(fit_ml_online)[names(fit_ml_online) == "se_ml_online"] <- "se"

fit_nls_online <- avg.pred.diff.nls_online$fit
se_nls_online <- avg.pred.diff.nls_online$se.fit

model_type <- c("NLS", "NLS", "NLS")

fit_nls_online <- cbind(fit_nls_online, se_nls_online, interview, model_type)
names(fit_nls_online)[names(fit_nls_online) == "se_nls_online"] <- "se"

fit_mle_online <- avg.pred.diff.mle_online$fit
se_mle_online <- avg.pred.diff.mle_online$se.fit

model_type <- c("ML with Floor/Ceiling", "ML with Floor/Ceiling", "ML with Floor/Ceiling")

fit_mle_online <- cbind(fit_mle_online, se_mle_online, interview, model_type)
names(fit_mle_online)[names(fit_mle_online) == "se_mle_online"] <- "se"

predplots_diff_online.df <- rbind(fit_ml_online, fit_nls_online, fit_mle_online)


ggplot(data = predplots_diff_online.df, aes(x=interview, y=fit, ymin=lwr, ymax=upr, colour=model_type, group=model_type)) + 
  expand_limits(y=c(0,1)) + 
  geom_point(size=8, stat="identity", position=position_dodge(0.4)) +
  geom_errorbar(aes(), width=.1, position=position_dodge(0.4)) +
  scale_colour_manual(values=c("#2c7bb6",  "#5e3c99", "#d7191c"), name="Model Type") + 
  scale_x_discrete(name="\nExperimental Condition",
                   limits=c("White Interviewer", "Latinx Interviewer", "Difference"))+  
  ylab("Estimated Proportion") +
  geom_text(aes(label=format(fit, digits=2)), fontface="bold", colour="white", size=3, position=position_dodge(0.4))+
  geom_hline(yintercept=0, linetype="dashed", color = "red") +
  theme_light(base_size=12) +
  theme(legend.position = c(.9, .75), legend.justification = c(.8, .2))


##Table A16. Online Survey, by Respondent Gender
reg_result_online_women <- lm(list ~ treat*int_hisp + age + college_grad + int_fem, data = women_online)
summary(reg_result_online_women)

reg_result_online_men <- lm(list ~ treat*int_hisp + age + college_grad + int_fem, data = men_online)
summary(reg_result_online_men)

stargazer(reg_result_online_men, reg_result_online_women)

reg_result_online_genderinteract <- lm(list ~ treat*int_hisp*female + age + college_grad + int_fem, data = online)
summary(reg_result_online_genderinteract)

vars.order <- c("treat", "int_hisp", "treat:int_hisp", "age", "college_grad", "int_fem", "female", "treat:female", "int_hisp:female", "treat:int_hisp:female")
stargazer(reg_result_online_men, reg_result_online_women, reg_result_online_genderinteract,
          order=paste0("^", vars.order , "$"),
          column.labels = c("\\thead{Male Participants}", "\\thead{Female Participants}", "\\thead{Triple Interaction}"),
          #type="latex", out="Tables/yougov_summary.tex", float = FALSE, #style="apsr",
          covariate.labels = c("Immigration Condition", "Latinx Interviewer", "Immigration Condition*Latinx Interviewer", "Age", "College Graduate", "Female Interviewer", "Female Participant", "Immigration Condition*Female Participant", "Latinx Interviewer*Female Participant", "Immigration Condition*Latinx Interviewr*Female Participant") 
)




reg_result_online_women = emmeans(reg_result_online_women, specs=revpairwise ~ treat:int_hisp)
reg_result_online_women$emmeans
reg_result_online_women$contrasts %>% summary(infer=TRUE) %>% as.data.frame()

reg_result_online_men = emmeans(reg_result_online_men, specs=revpairwise ~ treat:int_hisp)
reg_result_online_men$emmeans
reg_result_online_men$contrasts %>% summary(infer=TRUE) %>% as.data.frame()



##Table A17. Online Study: Linear Model with Interviewer Gender Interaction
reg_result_online_intfem <- lm(list ~ treat*int_fem + age + college_grad + female + int_hisp, data = online)
summary(reg_result_online_intfem)
stargazer(reg_result_online_intfem)

##Table A18. Difference-in-Differences Results by Interviewer Gender: Online Study
emm_online_intfem = emmeans(reg_result_online_intfem, specs=revpairwise ~ treat:int_fem)
emm_online_intfem$emmeans
emm_online_intfem$contrasts %>% summary(infer=TRUE) %>% as.data.frame()



##Table A19. Direct Immigration Attitudes: Linear Regression Model
##Face-to-Face
##In general, illegal immigrants are more prone to violence
reg_explicit1_inperson <- lm(Q11_1 ~ int_hisp + age + college_grad + female + int_fem, data = inperson)
summary(reg_explicit1_inperson)

##Given the current immigration situation, denying some basic constitutional rights is justified
reg_explicit3_inperson <- lm(Q11_3 ~ int_hisp + age + college_grad + female + int_fem, data = inperson)
summary(reg_explicit3_inperson)

##Online
##In general, illegal immigrants are more prone to violence
reg_explicit1_online <- lm(Q18_1 ~ int_hisp + age + college_grad + female + int_fem, data = online)
summary(reg_explicit1_online)

##Given the current immigration situation, denying some basic constitutional rights is justified
reg_explicit3_online <- lm(Q18_3 ~ int_hisp + age + college_grad + female + int_fem, data = online)
summary(reg_explicit3_online)


stargazer(reg_explicit1_inperson, reg_explicit1_online, reg_explicit3_inperson, reg_explicit3_online)


##Other Results
cohens_d(Q11_1 ~ int_hisp, data=inperson)
cohens_d(Q11_3 ~ int_hisp, data=inperson)

cohens_d(Q18_1 ~ int_hisp, data=subset(online, online$int_hisp!="NA"))
cohens_d(Q18_3 ~ int_hisp, data=subset(online, online$int_hisp!="NA"))

describe(inperson$Q11_1)
describe(inperson$Q11_3)
describe(online$Q18_1)
describe(online$Q18_3)

